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Reversible gelation and dynamical arrest of dipolar colloids 
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Abstract. - We use molecular dynamics simulations of a simple model to show that dispersions of 
slightly elongated colloidal particles with long-range dipolar interactions, like ferrofluids, can form 
a physical (reversible) gel at low volume fractions. On cooling, the particles first self-assemble 
into a transient percolating network of cross-linked chains, which, at much lower temperatures, 
then undergoes a kinetic transition to a dynamically arrested state with broken ergodicity. This 
transition from a transient to a frozen gel is characterised by dynamical signatures reminiscent of 
jamming in much denser dispersions. 
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Colloidal dispersions have been observed to undergo 
a glass transition to a disordered solid upon increasing 
their volume fraction <j> [1], or to undergo clustering and 
then gelation at low cf> upon increasing the mutual attrac- 
tion between particles [2]. Reversible gelation of colloidal 
dispersions is characterised by long-lived physical (rather 
than irreversible chemical) bonding between particles [3], 
which leads to the formation of sample-spanning networks 
capable of sustaining weak stresses at low volume frac- 
tions. In recent experimental and theoretical investiga- 
tions the gelation process is almost invariably associated 
with short-range attractive interactions between particles, 
induced, e.g., by the polymer-driven depletion mechanism, 
which lead to low coordination clustering or jamming of 
preformed clusters [2,4-7]. 

In this Letter we explore an alternative gelation mecha- 
nism driven by long-range anisotropic interactions result- 
ing from electric or magnetic dipoles as in the techno- 
logically important ferrofluids. Spherical particles with a 
point dipole are known to form a "string phase" at low vol- 
ume fraction, comprising long transient chains of particles 
concatenated in energetically favourable head-to-tail con- 
figurations [8] . Individual chains percolate in such a phase, 
i.e., the average chain length diverges and in simulations 
the particles become connected to their own periodic im- 
ages through the boundary conditions on the simulation 
cell. However, the chains do not appear to be sufficiently 
interconnected to form a genuine gel-like network. Here, 
we show that network formation may be achieved at low 4> 
by a modest elongation of the particles and their dipoles 



to make short dipolar dumbbells. 

Chain formation in dipolar colloids is a natural conse- 
quence of the anisotropy of the long-ranged dipole-dipole 
interactions. Colloidal chaining can also be achieved us- 
ing short-ranged attraction, by limiting each particle to 
a maximum of two (reversible) bonds at specified sticky 
patches on their surface. To form networks of such parti- 
cles, it is necessary to introduce a second type of particle 
with higher valency to act as chain junctions [7]. Alterna- 
tively, a low coordination number can be enforced by the 
introduction of three-body forces that discourage small 
angles between the bonds to any given particle [6]. In 
contrast, junctions arise with decreasing temperature in 
our one-component dipolar dumbbell fluid simply because 
their statistical weight increases relative to that of discon- 
nected chains. 

The nature [9] and location of any first-order fluid- 
fluid phase transition in dipolar spheres is still not com- 
pletely resolved. The critical point of the Stockmayer fluid 
(Lennard- Jones particles with a point dipole) has been 
shown to drop rapidly in temperature as the isotropic 
attraction is switched off, leaving soft dipolar spheres 
[10]. In contrast, Monte Carlo simulations of dipolar hard 
spheres indicate as many as three disordered phases at low 
temperature [11]. Less controversially, it has been shown 
that hard spherocylinders [12] or dumbbells [13] with an 
axial point dipole exhibit a gas-liquid transition within a 
limited range of particle elongations. However, the present 
simulations of soft-core dumbbells with an extended dipole 
showed no sign of such a transition within the range of 
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Fig. 1: Simulation snapshots at (a) T* = 0.205 and (b) T* = 
0.041. Dark, intermediate and light shading indicate particles 
with, respectively, one, two, and three or more neighbours. For 
simplicity, the slightly elongated dumbbells are represented by 
spheres. 



temperatures studied, as evidenced by the small wavevec- 
tor behaviour of the static structure factor discussed later 
in the Letter. 

Our dumbbells are made up of two interpenetrating col- 
loidal soft spheres carrying opposite charges ±q at their 
centres, separated by a fixed distance d. The resulting 
dipole moment is p = qd and the point dipole limit is 
recovered as d — > and q — > oo at fixed p. The inter- 
action potential between sites on different dumbbells is 
the sum of a soft-sphere repulsion c/r 12 , and a Coulomb 
interaction ±q 2 / Aire^er, where r is the site-site distance, 
while eo and e are the permittivity of free space and the 
relative permittivity of the medium, respectively. The 
finite extension of the dipole enhances string formation 
[14], while the dumbbell geometry facilitates branching of 
chains, which can thus interconnect into a genuine three- 
dimensional space-spanning network when the tempera- 
ture T is lowered at fixed <j). A convenient length scale is 
the distance a between the centres of two dumbbells in the 
head-to-tail configuration of lowest total potential energy, 
while the natural energy scale is the dipole-dipole energy 
u = pt 2 / 4ireQea 3 . In terms of these quantities, the soft- 
core repulsion coefficient is fixed at c = 0.0208uct 12 . We 
shall use the reduced units T* — ksT/u and p* = pa 3 
(cf> = irp*/6), where ks is Boltzmann's constant and 
p = N/V is the number of particles per unit volume, while 
the reduced time is t* = t/to with to = \J ma 1 ju and m 
the particle mass. For typical colloidal particles in water, 
a ~ 10 2 nm, m ~ 10~ 18 kg, and q ~ 10 2 proton charges so 
that for d/a = 0.217, as used in the present simulations, 
u/ks — 10 5 K and to ~ 10 _7 s, while room temperature 
corresponds to T* ~ 0.01. 

The constant-temperature molecular dynamics (MD) 
simulations were carried out on periodic samples of N = 
1000 dipolar dumbbells, using the Berendsen thermo- 
stat [15] and the standard Verlet leap-frog algorithm [16], 
implemented in the Gromacs package [17]. Long-range 
Coulomb interactions are dealt with by the particle-mesh 
Ewald method. The standard time step was chosen to 




rla 



Fig. 2: Radial distribution functions, g(r) at four reduced tem- 
peratures T*. Inset: probability distribution P(n) of the num- 
ber of neighbours n per particle. Lines join points at integer 
values of n for clarity. 



be At* = 0.0008 and most simulation runs extended over 
roughly 10 6 time steps after initial equilibration, with a 
few runs of up to ten times this length at the lowest tem- 
peratures. The Newtonian dynamics used here neglect the 
stochastic (Brownian) and hydrodynamic forces induced 
by the solvent. No single technique without explicit sol- 
vent fully captures these effects, and Newtonian dynam- 
ics remains appealing because of its efficiency and lack of 
ambiguity. It has also been shown, at least for denser sys- 
tems, that long-time dynamical evolutions are not affected 
by the precise short-time dynamics [18, 19]. 

All simulations were carried out along the isochore 
(f> = 0.0745 and T* was gradually lowered from T* ~ 0.4 
down to T* ~ 0.02. Typical snapshots of MD-generated 
configurations are shown in fig. 1 for T* — 0.205 and 
T* = 0.041. The snapshots highlight singly-connected 
(end-of-chain) , doubly-connected (mid-chain), and triply- 
or more highly-connected (junction) particles. We deter- 
mine the connectivity of two given particles by the simple 
geometric criterion that a pair of opposite charges — one on 
each dumbbell — lie closer than r cut = a. The higher tem- 
perature snapshot shows long chains which are hardly con- 
nected, while the lower temperature configuration is typ- 
ical of a percolating network. The percolation threshold, 
where 50% of configurations contain a system-spanning 
cluster, is found to occur at T* — 0.177. At this tem- 
perature, the network is highly transient and the chains 
are only loosely connected. Slightly below T p , the fraction 
of percolating configurations rapidly rises to unity as the 
largest cluster becomes more ramified. Nevertheless, the 
lifetime of individual connections between chains exceeds 
typical simulation time scales (i.e., millions of time steps) 
only at much lower temperatures. 

The temperature dependence of the pair distribution 
function g(r) of the particle centres is illustrated in fig. 2. 
At the highest temperature, a single peak appears at 
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r — <r, signalling the formation of dimers. As T* drops, 
the amplitude of the main peak rises sharply and well- 
separated secondary peaks develop at roughly multiple 
values of a as would be the case for a one-dimensional 
system, pointing to chain formation. The overall form 
of g(r) differs considerably from that of radial distribu- 
tion functions for dense, highly-coordinated fluids. The 
inset shows the probability distribution of the number n 
of neighbours per particle. Above T* most particles have 
a single or no neighbour, while below T* almost all par- 
ticles are doubly coordinated and a significant fraction of 
particles have three neighbours, indicating the existence of 
branched chains, which are crucial for network formation. 

The Fourier transform of g(r) yields the structure fac- 
tor S(q) which is experimentally accessible by X-ray or 
small-angle neutron diffraction. Except at the highest 
temperatures S(q) exhibits a pronounced low-q peak, in- 
dicative of clustering, similar to that observed for other 
models of gelation [6] . Extrapolation of the MD data to 
q = provides an estimate of the isothermal compress- 
ibility xt via the standard fluctuation formula [20]. The 
amplitude of the central peak implies that the gel is signif- 
icantly more compressible than the corresponding gas of 
non-interacting particles. S(q — 0) appears to grow slowly 
as T drops, but remains finite, precluding the possibility 
of spinodal instability associated with phase separation. 
The pressure is readily calculated from the virial theorem 
and turns out to drop rapidly with T due to clustering, as 
expected. 

Another important static characteristic is the relative 
dielectric permittivity e, which may be estimated from the 
fluctuations of the total dipole moment M of the sample 
using Kirkwood's relation [20] adapted to metallic bound- 
ary conditions: 
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e exhibits a striking non-monotonic behaviour as a func- 
tion of temperature. It first increases as T drops, due to 
the 1/T factor in eq. (1) and to enhanced dipolar fluctu- 
ations associated with cluster formation, e then reaches 
a maximum e ~ 10 at T* ~ 0.2, before dropping sharply 
to values close to the ideal gas value e = 1 at the lowest 
temperatures. Such low values show that dipolar fluctu- 
ations become strongly suppressed at very low tempera- 
tures, indicating incipient dynamical arrest. At such tem- 
peratures, the topology of connections between the col- 
loidal chains "freezes" into a quasi-permanent arrange- 
ment, thereby restricting the motion of particles to the 
scope permitted by deformations of the chains. 

The expected dynamical slowing down at low tempera- 
tures is best characterised by time-dependent correlation 
functions and transport coefficients. We have mostly fo- 
cused on single-particle motion which is sampled with N 
times better statistics than collective dynamics in the sim- 
ulations. Examples of the normalised velocity autocorre- 
lation function Z(t) = (v(t) ■ v(0)) / (\v\ 2 ) , where v(t) is 
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Fig. 3: Velocity autocorrelation function Z(t*). The arrow 
intersects the curves in order of increasing temperature. 



the velocity of a given particle at time t, are shown in 
fig. 3 for several temperatures. As expected, Z(t*) is seen 
to relax more and more rapidly as T* decreases. However, 
just below the percolation threshold at T* = 0.177, Z(t*) 
begins to exhibit a striking oscillatory behaviour. Oscilla- 
tions in Z(t) normally involve periods of negative values, 
which signal the backscattering of particles in a dense fluid 
as the particles, on average, reverse their velocity. In con- 
trast, our low-density gel exhibits a Z(t) that oscillates 
as it goes to zero without changing sign, except at the 
very lowest T*. Two frequencies can be discerned. The 
shorter period (0.3 to 0.5io) can be attributed to the one- 
dimensional rattling of particles within the chains. This 
assignment is compatible with the Einstein frequency de- 
rived from a Taylor expansion to the quadratic term [21] of 
the potential energy of a dumbbell moving in one dimen- 
sion between two fixed dumbbells with aligned dipoles. 
This rapid oscillation is superimposed upon a slower move- 
ment of period 1.8to corresponding to concerted movement 
of a chain segment [22]. 

Figure 4(a) shows the mean-square displacement (MSD) 
of a particle from its original position, (\r(t) — r(0)| 2 ) ver- 
sus time, on a double- logarithmic plot. At the higher tem- 
peratures the MSD goes over directly from the initial bal- 
listic regime (~ t 2 ) to the linear Einstein regime at longer 
times [22]. Below T* a sub-diffusive regime, characterised 
by a clear inflection point sets in at intermediate times. 
The sub-diffusive regime extends over a rapidly increas- 
ing time interval as T* drops, and at the lowest T* the 
normal diffusive regime is only just reached within the 
time-scale of our simulation. In other words, the system 
becomes dynamically arrested, at least over the accessible 
time scale, which covers several decades. The particle dif- 
fusion constant D* is given by the slope of the MSD in 
the linear regime or by the time integral of Z(t) [20] and 
drops rapidly with T*, as illustrated in fig. 4(b). At low 
temperatures, where diffusion is an activated process, D* 
follows an Arrhenius law. 
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Fig. 4: (a) Reduced mean-square displacement (MSD) as a 
function of time. Dashed lines indicate log-log slopes of 2 
(ballistic regime) and 1 (diffusive). The arrow intersects the 
curves in order of increasing temperature, (b) Arrhenius plot 
of the reduced diffusion constant D* (circles) and reduced rate 
constant kl for bond breaking (squares). 



Another instructive diagnostic, also plotted in fig. 4(b), 
is provided by the rate of bond breaking. Denning bij (t) to 
be unity if dumbbells i and j are neighbours at time t and 
zero otherwise, we may construct a bond correlation func- 
tion B(t) — (b(t)b(0)), which exhibits a regime of expo- 
nential decay over a wide range of temperatures, allowing 
a reduced rate constant k£ to be defined. As seen for D*, 
fc^ decreases rapidly upon cooling from high temperatures, 
and enters a regime of Arrhenius temperature dependence 
once the network is properly established. However, at very 
low temperatures another change is seen, and k£ is higher 
than an extrapolation would predict. The relative ease of 
breaking bonds here can be understood by noting that a 
dipolar particle can only make two energetically optimal 
bonds (i.e., with dipoles aligned head-to-tail). The in- 
creasing average coordination number at low temperature 
necessarily means that many particles are forming bonds 
with frustrated geometries, which require less energy to 
break. 

The local motion of individual particles is characterised 
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Fig. 5: Self-intermediate scattering function F s (q,t) at qo = 
1.12 as a function of time scaled by the relaxation time r s at 
each temperature. Inset: F s (q,t) at T* — 0.020 as a function 
of reduced time t* . 



by the the self-intermediate scattering function of incoher- 
ent neutron scattering experiments [20] 
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where q is a wavevector compatible with the periodic 
boundary condition in the simulations. A g-dependent 
relaxation time r s (g) may be defined by the time- integral 
of (2). At high temperatures, a ballistic regime is ap- 
proached and F s (q,t) shows Gaussian decay [20]. This is 
illustrated in fig. 5, where F s (q, t) is plotted versus t/r s (q) 
for qo = 1.12. For T < T*, the correlation function devi- 
ates increasingly from the Gaussian limit, but never con- 
forms to the exponential decay that would be expected in 
a dense liquid. At sufficiently low temperature, a stretched 
long-time tail develops. At the lowest temperature investi- 
gated, T* = 0.020 (inset), F 8 (q,t) exhibits strong dynam- 
ical slowing down, with an initial relaxation, followed 
by a plateau (which defines a "non-ergodicity" parame- 
ter f q ~ 0.895) and a final a-relaxation for t* > 170, 
using the standard terminology of the kinetic glass tran- 
sition [20,23]. Qualitatively similar behaviour is observed 
for the collective intermediate scattering function (or den- 
sity autocorrelation function) F(q,t), although the latter 
is subject to much larger statistical noise, and requires 
averages to be taken over several replicas, i.e., different 
initial conditions of the system. 

The dependence of the self- intermediate scattering func- 
tion on the wavevector q is illustrated in fig. 6 at T* = 
0.164 and T* — 0.020 (inset). The former temperature 
lies just below T*, where a regime of exponential decay is 
rapidly reached at all wavevectors, indicating that the pre- 
liminary formation of the network is not associated with 
any dynamical arrest. The initial curvature of the relax- 
ation at large q is due to the very brief ballistic travel of 
particles before they reach their neighbour. This region 
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Fig. 6: Dependence of the self-intermediate scattering function 
F s (q,t) on the wavevector q at T=0.164 as a function of time 
scaled by the relaxation time r s at each wavevector. Inset: q- 
dependence of F a (q,t) at T* — 0.020 as a function of reduced 
time t* . Arrows associate each curve with either the left-hand 
or the right-hand vertical scale. 

is emphasised for large q because the time axis has been 
scaled by r s (q), which decreases rapidly with q. In con- 
trast, at the lower temperature (inset of fig. 6), F s (q,t) 
encounters a g-dependent plateau, typical of glassy dy- 
namics and long-lived dynamical arrest. 

The data shown in the inset of fig. 6 are semi- 
quantitatively reproducible with different initial condi- 
tions, hinting that full dynamical arrest may only occur 
in the zero-temperature limit in low-density gels, contrary 
to the case of high-density structural glasses which ex- 
hibit clear-cut non-Arrhenius behaviour at finite temper- 
atures. This qualitative difference in behaviour may be at- 
tributed to the fact that particles are not completely caged 
by neighbours in the open network, and that "floppy" mo- 
tions of the latter are not strongly inhibited at low packing 
fractions. 

In summary our simulations predict that dilute disper- 
sions of dumbbell-shaped, dipolar colloidal particles un- 
dergo a two-step gelation transition, with a reversible, 
ergodic network gel formed at a percolation threshold, 
characterised by relatively short-lived bonds, as in "liv- 
ing" polymers [24], followed by a kinetic transition to a 
"frozen", ergodicity-breaking gel at lower temperatures. 
The predicted dynamical behaviour is not unlike that ob- 
served in very different colloidal systems involving short- 
range attractive forces between particles [2,4-7], thus hint- 
ing at a certain universality of structural arrest for col- 
loidal gels. 

The model put forward in this letter could be mimicked 
experimentally using suspensions of sterically stabilised, 
oppositely charged colloidal particles assembled in pairs 
by the emulsion technique developed by Pine and collabo- 
rators [25,26]. The wide variation of T* could be achieved 



by tuning the particle dipole moment /i, e.g., by vary- 
ing the surface charge. Future work will investigate the 
dependence of the present gelation mechanism on the col- 
loidal volume fraction, which is a more convenient control 
parameter in experiments. 
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